% ad.m

function [Gz,Gzp,Gnp,Fz,Fzp,Fnp, ... 
    GFzz,GFzpz,GFnpnp,GFznp,GFzpnp,GFzpzp, ... 
    GFzpzpzp,GFzpzpnp,GFzpzpz,GFzpzz,GFzpznp,GFzpnpnp,GFzzz,GFzznp,GFznpnp,GFnpnpnp] ... 
    =adm(F,G,z,zp,np,varargin)
% Syntax
% =======
% 
% [Gz,Gzp,Gnp,Fz,Fzp,Fnp, 
%    GFzz,GFzpz,GFnpnp,GFznp,GFzpnp,GFzpzp, 
%    GFzpzpzp,GFzpzpnp,GFzpzpz,GFzpzz,GFzpznp,GFzpnpnp,GFzzz,GFzznp,GFznpnp,GFnpnpnp] 
%    =adm(F,G,z,zp,np,...)
% 
% Options
% ========
% 
% - 'cache=': re-use symbolic derivatives from the previous call if there
% are no changes to model equations or options. Default: true. 
% - 'parallel=': use the MATLAB Parallel Computing Toolbox. Default: false.
% - 'order=': order of differentiation. Default: 2. 
% 
% Description
% ============
% 
% Computes arrays symbolic derivatives as defined in Johnston, King, and 
% Lie (2014) using the MATLAB Symbolic Math Toolbox. 

persistent oF oG oz ozp onp ovarargin ...
    oFz oFzp oFnp oGz oGzp oGnp ...
    oGFzz oGFzpz oGFnpnp oGFznp oGFzpnp oGFzpzp ...
    oGFzpzpzp oGFzpzpnp oGFzpzpz oGFzpzz oGFzpznp oGFzpnpnp oGFzzz oGFzznp oGFznpnp oGFnpnpnp

[cache,varargin] = parseOption(varargin,'cache',true) ;
[order,varargin] = parseOption(varargin,'order',2) ;
[parallel,varargin] = parseOption(varargin,'parallel',false) ;

nz = length(z) ;
nn = length(np) ;
nf = length(F) ;
ng = length(G) ;

assert(nf+ng == nz, 'Number of equations must equal number of variables.') ;
assert(nz == size(zp,2), 'Variable vectors are constructed inconsistently.') ;

if cache
    if not(isempty(oF)) && not(isempty(oG))
        if all(size(F)==size(oF)) && all(size(G)==size(oG))
            try 
                crit = norm([eval(F-oF);eval(G-oG)]) ;
            catch
                crit = Inf ;
            end
        else crit = Inf ;
        end
        if numel(varargin)==numel(ovarargin)
            veq = true ;
            for ii = 1:numel(varargin)
                switch class(varargin{ii})
                    case 'char'
                        if ~strcmpi(varargin{ii},ovarargin{ii}), veq = false ; end
                    otherwise
                        if varargin{ii}~=ovarargin{ii}, veq = false; end
                end
            end
        else veq = false ;
        end
        if not(veq), crit = Inf ; end
        if crit<sqrt(eps) && strcmp(char(z),char(oz)) && strcmp(char(zp),char(ozp)) && strcmp(char(np),char(onp))
            fprintf(1,'No change to input arguments detected, using cached derivatives...\n') ;
            Fz = oFz ;
            Fzp = oFzp ;
            Fnp = oFnp ;
            Gz = oGz ;
            Gzp = oGzp ;
            Gnp = oGnp ;
            GFzz = oGFzz ;
            GFzpz = oGFzpz ;
            GFnpnp = oGFnpnp ;
            GFznp = oGFznp ;
            GFzpnp = oGFzpnp ;
            GFzpzp = oGFzpzp ;
            GFzpzpzp = oGFzpzpzp ;
            GFzpzpnp = oGFzpzpnp ;
            GFzpzpz = oGFzpzpz ;
            GFzpzz = oGFzpzz ;
            GFzpznp = oGFzpznp ;
            GFzpnpnp = oGFzpnpnp ;
            GFzzz = oGFzzz ;
            GFzznp = oGFzznp ;
            GFznpnp = oGFznpnp ;
            GFnpnpnp = oGFnpnpnp ;
            return
        end % if crit<sqrt(eps)
    end
end

if order>1
    GFzz = sym(zeros(nz,nz^2)) ;
    GFzpz = sym(zeros(nz,nz^2)) ;
    GFzpzp = sym(zeros(nz,nz^2)) ;
    GFnpnp = sym(zeros(nz,nn^2)) ;
    GFznp = sym(zeros(nz,nn*nz)) ;
    GFzpnp = sym(zeros(nz,nn*nz)) ;
    if order>2
        GFzpzpzp = sym(zeros(nz,nz^3)) ;
        GFzpzpz = sym(zeros(nz,nz^3)) ;
        GFzpzpnp = sym(zeros(nz,nz^2*nn)) ;
        GFzpzz = sym(zeros(nz,nz^3)) ;
        GFzpznp = sym(zeros(nz,nz^2*nn)) ;
        GFzpnpnp = sym(zeros(nz,nn^2*nz)) ;
        GFzzz = sym(zeros(nz,nz^3)) ;
        GFzznp = sym(zeros(nz,nz^2*nn)) ;
        GFznpnp = sym(zeros(nz,nn^2*nz)) ;
        GFnpnpnp = sym(zeros(nz,nn^3)) ;
    end
end

%% Differentiation
fprintf('\nConstructing analytical derivatives... ') ;

% First order:
Fz = jacobian(F,z) ;
Fzp = jacobian(F,zp) ;
Fnp = jacobian(F,np) ;
Gz = jacobian(G,z) ;
Gzp = jacobian(G,zp) ;
Gnp = jacobian(G,np) ;

if order>1
    mndiff = @(X,v) transpose( vec( transpose( jacobian(X,v) ) ) ) ;
    
    GFz = [Gz; Fz] ;
    GFzp = [Gzp; Fzp] ;
    GFnp = [Gnp; Fnp] ;
    
    % Second order:
    if parallel
        parfor h = 1:nz
            GFzz(h,:) = mndiff(GFz(h,:),z) ;
            GFzpz(h,:) = mndiff(GFzp(h,:),z) ;
            GFzpzp(h,:) = mndiff(GFzp(h,:),zp) ;
            GFnpnp(h,:) = mndiff(GFnp(h,:),np) ;
            GFznp(h,:) = mndiff(GFz(h,:),np) ;
            GFzpnp(h,:) = mndiff(GFzp(h,:),np) ;
        end
    else
        for h = 1:nz
            GFzz(h,:) = mndiff(GFz(h,:),z) ;
            GFzpz(h,:) = mndiff(GFzp(h,:),z) ;
            GFzpzp(h,:) = mndiff(GFzp(h,:),zp) ;
            GFnpnp(h,:) = mndiff(GFnp(h,:),np) ;
            GFznp(h,:) = mndiff(GFz(h,:),np) ;
            GFzpnp(h,:) = mndiff(GFzp(h,:),np) ;
        end
    end
    if order>2
        % Third order
        if parallel
            parfor h = 1:nz
                GFzpzpzp(h,:) = mndiff(GFzpzp(h,:),zp) ;
                GFzpzpz(h,:) = mndiff(GFzpzp(h,:),z) ;
                GFzpzpnp(h,:) = mndiff(GFzpzp(h,:),np) ;
                GFzpzz(h,:) = mndiff(GFzpz(h,:),z) ;
                GFzpznp(h,:) = mndiff(GFzpz(h,:),np) ;
                GFzpnpnp(h,:) = mndiff(GFzpnp(h,:),np) ;
                GFzzz(h,:) = mndiff(GFzz(h,:),z) ;
                GFzznp(h,:) = mndiff(GFzz(h,:),np) ;
                GFznpnp(h,:) = mndiff(GFznp(h,:),np) ;
                GFnpnpnp(h,:) = mndiff(GFnpnp(h,:),np) ;
            end
            
        else
            for h = 1:nz
                GFzpzpzp(h,:) = mndiff(GFzpzp(h,:),zp) ;
                GFzpzpz(h,:) = mndiff(GFzpzp(h,:),z) ;
                GFzpzpnp(h,:) = mndiff(GFzpzp(h,:),np) ;
                GFzpzz(h,:) = mndiff(GFzpz(h,:),z) ;
                GFzpznp(h,:) = mndiff(GFzpz(h,:),np) ;
                GFzpnpnp(h,:) = mndiff(GFzpnp(h,:),np) ;
                GFzzz(h,:) = mndiff(GFzz(h,:),z) ;
                GFzznp(h,:) = mndiff(GFzz(h,:),np) ;
                GFznpnp(h,:) = mndiff(GFznp(h,:),np) ;
                GFnpnpnp(h,:) = mndiff(GFnpnp(h,:),np) ;
            end
        end
    end
end

%% Handle empty output arguments
if order<3
    GFzpzpzp = sym(0) ;
    GFzpzpnp = sym(0) ;
    GFzpzpz = sym(0) ;
    GFzpzz = sym(0) ;
    GFzpznp = sym(0) ;
    GFzpnpnp = sym(0) ;
    GFzzz = sym(0) ;
    GFzznp = sym(0) ;
    GFznpnp = sym(0) ;
    GFnpnpnp = sym(0) ;
    if order<2
        GFzz = sym(0) ;
        GFzpz = sym(0) ;
        GFnpnp = sym(0) ;
        GFznp = sym(0) ;
        GFzpnp = sym(0) ;
        GFzpzp = sym(0) ;
    end
end

%% Handle caching, if appropriate
if cache
    oF = F ;
    oG = G ;
    oz = z ;
    ozp = zp ;
    onp = np ;
    ovarargin = varargin ;
    oFz = Fz ;
    oFzp = Fzp ;
    oFnp = Fnp ;
    oGz = Gz ;
    oGzp = Gzp ;
    oGnp = Gnp ;
    oGFzz = GFzz ;
    oGFzpz = GFzpz ;
    oGFnpnp = GFnpnp ;
    oGFznp = GFznp ;
    oGFzpnp = GFzpnp ;
    oGFzpzp = GFzpzp ;
    oGFzpzpzp = GFzpzpzp ;
    oGFzpzpnp = GFzpzpnp ;
    oGFzpzpz = GFzpzpz ;
    oGFzpzz = GFzpzz ;
    oGFzpznp = GFzpznp ;
    oGFzpnpnp = GFzpnpnp ;
    oGFzzz = GFzzz ;
    oGFzznp = GFzznp ;
    oGFznpnp = GFznpnp ;
    oGFnpnpnp = GFnpnpnp ;
end

fprintf('%g\n', toc) ;

end






